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Abstract 

Management  of  the  water  and  heat  ejected  as  byproducts  in  an  operating  PEM  fuel  cell  stack  are  crucial  factors  in  their  optimal  design  and  safe 
operations.  Models  currently  available  for  a  PEM  fuel  cell  are  based  on  either  empirical  or  3-D  computational  fluid  dynamics  (CFD).  Both  models 
do  not  fully  meet  the  need  to  represent  physical  behavior  of  a  stack  because  of  either  their  simplicity  or  complexity.  We  propose  a  highly  dynamic 
PEM  fuel  cell  stack  model,  taking  into  account  the  most  influential  property  of  temperature  affecting  performance  and  dynamics.  Simulations  have 
been  conducted  to  analyze  start-up  behaviors  and  the  performance  of  the  stack  in  conjunction  with  the  cells.  Our  analyses  demonstrate  static  and 
dynamic  behaviors  of  a  stack.  Major  results  presented  are  as  follows:  (1)  operating  dependent  temperature  gradient  across  through-plane  direction 
of  the  fuel  cell  stack,  (2)  endplate  effects  on  the  temperature  profile  during  start-up  process,  (3)  temperature  profile  influences  on  the  output  voltage 
of  individual  cells  and  the  stack,  (4)  temperature  influence  on  the  water  content  in  membranes  of  different  cells,  and  (5)  cathode  inlet  relative 
humidity  influence  on  the  temperature  profile  of  the  stack. 

©  2005  Elsevier  B.V.  All  rights  reserved. 
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1.  Introduction 

The  PEM  fuel  cell  is  a  potential  candidate  for  use  as  an 
alternative  power  source  in  future  vehicle  and  power  condition¬ 
ing  applications.  The  electric  power,  as  a  load,  continuously 
varies  with  a  function  of  time  in  the  systems.  Accordingly,  the 
flow  rates  of  fuels  and  reactants  should  be  properly  supplied  to 
meet  the  dynamics  of  load  requirement.  Hence,  temperature  and 
humidity  for  the  stack  should  be  precisely  maintained  to  secure 
an  efficient,  safe,  and  durable  operation  of  the  PEM  fuel  cell 
system. 

The  complexity  of  these  interrelated  physical  phenomena 
impedes  design  and  system  engineers  from  optimizing  design 
parameters  and  criterion  that  lead  to  high  costs  and  engineering 
time.  Utilization  of  computer  models  have  been  proven  to  be  the 
most  effective  tool  requiring  a  comprehensive  physical  model 
reflecting  real  behaviors  of  the  fuel  cell  stack. 
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Currently,  models  can  be  classified  into  two  groups:  a  simple 
composition  of  passive  electric  circuit  representing  the  double 
layer  effects  with  a  voltage  source  using  I-V  empirical  curve 
[1,2],  and  a  relatively  precise  cell  model  based  on  CFD  [3-5] 
whose  behavior  is  governed  by  equations  describing  chemi¬ 
cal  reactions,  fluid  and  thermodynamic,  and  electric  properties. 
These  CFD  models  take  enormous  computational  time  necessary 
for  a  calculation  of  grid  points  and  constrain  the  wide  use  of  the 
models.  Moreover,  associated  parameters  and  variables  required 
for  the  models  are  hard  to  characterize  and  measure  because  of 
the  thinness  of  layers  during  operations.  Thus,  these  models  are 
mostly  applied  to  investigate  parts  of  complex  domains  such  as 
flow  fields  of  a  single  cell  and  predict  the  cell  performance  with  a 
constant  temperature.  A  possible  loophole  has  been  proposed  by 
other  authors  with  lumped  models  emphasizing  dynamic  behav¬ 
ior  of  the  stack.  For  example,  Amphlett  et  al.  [6]  developed  a  first 
dynamic  model  with  empirical  thermal  values,  and  Gurski  et  al. 

[7]  expanded  it  to  consider  reactant  flow  and  coolant  control. 
Further  improvement  was  made  by  Muller  and  Stefanopoulou 

[8] ,  who  calculated  the  temperature  variation  of  a  stack.  How¬ 
ever,  the  simulation  results  do  not  incorporate  either  the  dynamic 
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2.  PEMFC  stack  modeling 

2.7.  Assumptions  and  domain 


Nomenclature 

A  area  (m2) 

C  mass  concentration  (kg  m- 3 ) 

Cp  specific  heat 

F  Faraday  number 

I  current  density  (A  m-2) 

/  thickness  (M) 

m  mass  (kg) 

M  mole  mass  (kg  mol-1) 

P  pressure  (partial  pressure)  (Pa) 

Rq\q  electrical  resistance  (£2) 

S  entropy  (J  mol- 1  K- 1 ) 

T  temperature  (K) 

Superscripts  and  subscripts 
cv  control  volume 

g  gas 

i,  j  index 

1  liquid 

an  anode 

ca  cathode 

sat  saturation 

sou  source 


or  transient  aspects  of  the  fuel  cell  system  in  operating  environ¬ 
ments. 

These  drawbacks  have  been  significantly  complemented  by 
using  simplified  thermodynamic  models  to  analyze  the  perfor¬ 
mance  of  the  stack,  by  Sundaresan  and  Moore  [9].  His  model 
based  on  layers  is  used  to  analyze  the  start-up  behavior  from 
a  sub-freezing  temperature.  However,  the  model  does  not  fully 
consider  several  factors:  (1)  flow  of  species  at  the  inlet  must 
be  the  same  as  that  at  the  outlet.  Thus,  no  fluid  dynamics  are 
considered  in  the  model.  (2)  Heat  source  terms  in  both  the 
catalysts  are  empirically  calculated  with  values  suggested  by 
Wohr  and  Peinecke’s  modeling  [10,1 1].  Accordingly,  the  anode 
source  term  is  presumed  as  a  relatively  large  value  that  should 
be  referred  to  be  around  zero  [12].  As  a  result,  the  model  does 
not  show  asymmetric  phenomena  of  performance  across  the 
stack. 

Wetton  et  al.  [13]  proposed  an  explicit  thermal  model  with 
the  coolant  channel  coupled  with  a  1-D  cell  model  [14],  which 
shows  an  outstanding  temperature  gradient  of  the  stack,  but  with 
no  dynamics  at  all. 

In  fact,  dehydration  of  membranes  dynamically  varies  with 
temperature,  which  strongly  influences  overall  performance 
of  cells  and  finally  the  resulting  stack.  Consequently,  a  high 
dynamic  model  for  a  stack  based  on  layers  of  a  single  cell  [15] 
has  been  developed  by  reflecting  the  dynamics  under  different 
operating  conditions  for  pressure,  humidity,  reactants  and  tem¬ 
perature.  Temperature  profiles  in  a  stack  with  associated  effects 
on  the  stack  performance  show  the  high  dependency  in  not  only 
static  but  also  dynamic  behavior. 


Assumptions  made  are  as  follows: 

1.  reactants  as  ideal  gases; 

2.  no  pressure  gradient  between  the  anodic  and  cathodic  side, 
which  results  in  no  convection,  but  only  diffusion  for  gas 
transport; 

3.  identical  inlet  conditions  of  each  cell  for  both  the  cathode 
and  anode  as  well  as  coolant  channel; 

4.  no  gas  pressure  drop  along  the  gas  channel; 

5.  linear  temperature  gradient,  linear  across  the  layers  in  the 
stack; 

6.  constant  thermal  conductivity  of  the  materials  in  a  fuel  cell; 

7.  no  contact  resistance; 

8.  negligible  anodic  over-potential; 

9.  no  current  density  gradient  across  the  cathode  catalyst  layer, 
which  implies  a  complete  reaction  of  reactants  at  the  cathode 
catalyst  layer  surface. 

10.  the  anodic  stoichiometrical  coefficient  is  1,  which  indicates 
a  complete  reaction  of  the  fuel  filled  in  the  fuel  cell  stack. 

Based  on  these  assumptions,  a  layered  based  model  for  a 
10-cells  stack  is  developed.  The  schematic  configuration  of  the 
simulation  domain  for  the  stack  is  shown  in  Fig.  1,  where  the 
cathode  sides  of  the  cells  are  located  on  the  left-hand  side,  while 
the  anodic  on  the  right-hand  side. 

2.2.  Single  cell 

2.2.1.  Governing  equations 

A  model  for  a  single  cell  has  been  developed  with  composi¬ 
tions  of  individual  layers  and  the  performance  detailed  analyzed 
[15].  Table  1  summarizes  the  sub-models  the  governing  the  equa¬ 
tions  describing  phenomenon  for  a  single  cell  [15]. 

2.2.2.  Thermal  model 

The  fuel  cell  is  sandwiched  between  two  bipolar  plates.  The 
individual  layer  of  a  single  cell  and  a  bipolar  plate  is  regarded  as  a 
control  volume,  whose  thermal-physical  properties  are  isotropic 
and  constant.  Then,  according  to  the  energy  conservation  equa¬ 
tion,  the  total  energy  changes  in  a  control  volume  are  equal 
to  the  sum  of  the  energy  exchange  at  boundaries  and  internal 

Table  1 


Equation  system  of  the  single  cell  model 


Sub-models 

Phenomenon 

Mathematical  equation 

Electrochemical  model 

Chemical  reaction 

Butler- Volmer 

Layered  thermal  model 

Temperature  variation 

Energy  conservation 

Membrane  water  balance 

Water  transport 

1-D  mass  conservation 

Catalyst  layer  proton  model 

Proton  concentration 

variation 

Empirical 

GDL  O2  diffusion  model 

Multi-component 

diffusion 

1-D  Stefan-Maxwell 

Gas  channel  mass  balance 

0-D  mass  conservation 
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Fig.  1.  Stack  schematic  configuration. 


energy  resources.  The  energy  exchanges  at  boundaries  occur 
by  three  factors:  (1)  mass  flow  into  each  volume;  (2)  conduc¬ 
tion  heat  transfer  through  the  cell;  (3)  convection  heat  transfer 
between  the  reactants  and  bipolar  plates  with  the  coolant.  Thus, 
the  thermal-dynamic  behavior  can  be  described  by  using  the 
following  energy  conservation  equation: 


E 


Op  Ci, mass  A cen /cv 


dTcv 
d  t 


i 


IP  in  A  cell  Op  (  Ein 


- v - 

Mass  flow  in 


T  GConv^-cell  T  Gcond^cell  T  Q 


sou 


■v* 


Convection  heat  transfer  Conduction  heat  transfer  Sources 


(1) 


Hence,  the  internal  energy  source  is  mainly  composed  of 
the  entropy  loss  and  the  chemical  energy  required  for  oxygen 
and  protons  to  overcome  the  barrier  of  the  over-potentials  in 
both  catalyst  layers  (Eq.  (2)).  In  addition,  other  heat  sources 
are  associated  with  ohmic  losses  by  a  transport  of  electrons  and 
protons  in  a  cell: 


.  T  AS 

2sou  —  iA cen  (  |-  f)  +  iAce\\ Rq\q 


In  fact,  the  change  in  entropy  due  to  the  electrochemical  reac¬ 
tion  (Eq.  (3)  and  (8))  in  both  of  the  catalysts  sides  predominantly 
influences  the  energy  sources  term,  which  can  be  clarified  by  the 
calculations  below. 

The  chemical  reaction  of  the  cathode  catalysts  can  be 
described  as, 


02(g)  +  2H(aq)  +  2e  ^  H20(i/g) 


Therefore,  the  entropy  change  of  the  cathode  reaction  is  equal 
to  the  sum  of  that  in  water,  oxygen  and  electron.  The  phase  of 
water  produced  can  be  either  vapor  or  liquid  which  makes  a 
difference  in  the  entropy  change  of  the  reaction.  The  value  of 
the  entropy  change  indirectly  measured  by  Ref.  [16]  amounts  to 
52  J  mol  - 1  K- 1 .  However,  this  value  is  quite  different  from  other 
authors  [3,17,18].  Therefore,  the  energy  generated  in  a  catalyst 
layer  for  three  cases  are  calculated  and  the  values  are  compared 
for  a  better  judgment. 

The  first  case  assumes  water  generated  as  vapor  phase.  Then, 
the  cathode  entropy  change  can  be  described  as  follows: 

9  1 

=  AlSca  =  s[H20(g)]  -  -402(g)]  -  24e“t] 

1 

=  198  -  -  x  205.03  -  2  x  65.29 
2 

=  -35.095  J  mol-1  K-1  (4) 


The  resulting  energy  by  the  entropy  change  is  equal  to: 
Q  =  A  ScaT  =  -35.0957  Jmor1  K"1 


The  second  case  assumes  water  as  liquid  phase.  The  corre¬ 
sponding  entropy  change  of  the  reaction  and  the  energy  result 
in  as  follows: 


nF 


dE 


o 


dr 


AScsl  =  4H20(i)] 


1 

-402(g)]  -  2s 


1 

=  69.91  -  -  x  205.03  -  2  x  65.29 
2 


=  —163  J  mol-1  K-1 


Q  =  A  ScaT  =  -163T  Jmol-1  K-1 
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The  last  case  assumes  that  the  water  generated  is  imme¬ 
diately  evaporated  at  the  catalysts.  The  corresponding  energy 
required  for  the  evaporation  can  be  expressed  as  a  product  of 
the  latent  heat  and  the  rate  of  mass  change  from  liquid  phase  to 
vapor: 

Q  —  hfgriifg  (7) 

A  comparison  with  the  experimental  value  [16]  indicates  that 
most  parts  of  water  generated  might  be  a  state  of  vapor.  Later 
in  Section  3,  simulations  have  been  conducted  for  the  first  two 
cases. 

In  contrast,  the  anodic  side  is  assumed  to  be  a  standard  elec¬ 
trode,  while  the  membrane  is  assumed  to  provide  a  water-like 
environment.  Thus,  the  anodic  entropy  change  becomes  almost 
zero. 

H2(g)  ^  ^H^q)  +  2e(pt)  (8) 


Fig.  2.  I-V  curve  for  different  cell  working  temperature. 


3.1.  Parameters 

The  parameters  used  for  models  are  summarized  as  follows 
(see  Table  2). 


2.3.  Stack 

The  model  for  a  stack  consists  of  models  for  endplates  assem¬ 
bly,  multiple  single  cells  and  coolant  channels  as  in  a  typically 
designed  stack.  The  stack  voltage  and  current  are  obtained  by 
series  and  parallel  connections  of  outputs  of  individual  cells, 
while  the  thermal  conditions  for  each  of  the  cells  and  coolant 
channels,  or  between  coolant  channels  and  endplate  assem¬ 
bly,  are  coupled  by  using  the  energy  conservation  equation 
(Eq.  (1)). 

3.  Simulation  results 

All  of  the  aforementioned  models  are  coded  with  Mat- 
lab/Simulink  and  C.  Simulations  with  different  operating  condi¬ 
tions  have  been  conveyed  to  investigate  the  static  and  dynamic 
behavior  of  a  single  as  well  as  a  stack,  which  include  the  typical 
polarization  curves  at  different  working  temperatures,  a  start-up 
behavior  and  transient  responses  of  the  stack  on  a  current  as  a 
step  load. 


3.2.  Polarization  curve 

Fig.  2  shows  the  temperature  dependent  I-V  characteristics 
from  333  to  363  K  with  a  step  of  10  K.  As  the  temperature 
increases,  the  water  removal  is  easier.  The  effects  are  consid¬ 
erably  high  at  the  range  of  the  higher  cell  current,  where  more 
water  is  produced.  This  result  is  comparable  with  an  analysis 
using  CFD  [19]. 

3.3.  Dynamic  behaviors 

Currently,  lack  of  experimental  data  on  the  water  phase  inside 
cells  impedes  any  statement  on  how  much  of  the  enthalpy  is  gen¬ 
erated  by  a  water  molecule,  which  includes  decisive  information 
on  the  state  of  water  generated.  Thus,  the  two  cases  of  either  the 
liquid  or  the  vapor  are  selected  and  summarized  as  follows: 

Case  1.  If  the  RH  of  the  cathode  inlet  is  low  and  at  the  same 
time  the  stoichiometrical  number  on  the  cathode  is  high,  water 
generated  by  the  chemical  reaction  in  the  catalysts  is  immedi¬ 
ately  removed  by  the  cathode  outlet  gas  in  the  form  of  vapor. 


Table  2 


Geometrical  parameters 


Thickness  (m) 

Heat  conductivity  (W  m  1  K  1 ) 

Heat  capacity  (J  kg  1  K  1 ) 

Density  (kg  m  3) 

GDL 

0.0004 

65 

840 

2000 

Catalyst  layer 

0.000065 

0.2 

770 

387 

Membrane  layer 

0.000183 

0.21 

1100 

1967 

Gas  channel 

0.001 

52 

935 

1400 

Plate 

0.001 

52 

935 

1400 

Coolant  channel 

0.001 

30 

935 

1400 

GDL  porosity 

0.5 

GDL  tortuosity 

3.725 

Bipolar  plate  contact  area 

percentage 

0.55 

Membrane  molecular  mass  (kg  mol  1 ) 

1.1 

Fuel  cell  area  (m2) 

0.0367 

Fuel  cell  active  area  (m2) 

0.03 
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Table  3 


Operating  condition  of  the  stack 


Gas  inlet  condition 

Pin 

(Pin 

RHin 

fo2  or  CO 

Cathode 

141000 

2.5 

0.2 

0.21 

Anode 

141000 

1 

0.001  or  0.5 

0 

This  state  can  be  regarded  as  a  vapor  phase,  where  no  liquid  is 
formed  in  the  stack.  One  of  possible  operating  conditions  is  listed 
in  Table  3  with  other  parameters  necessary  for  a  simulation: 

Case  2.  If  the  RH  at  the  anodic  inlet  is  high  and  stoichiometrical 
number  at  the  cathode  is  low,  then  water  generated  at  the  catalysts 
cannot  be  removed  immediately  in  the  form  of  vapor  in  the 
cathode  outlet  gas.  Thus,  assumptions  can  be  made  that  water  is 
mainly  generated  and  removed  in  the  liquid  phase.  These  cases 
are  listed  in  Table  4  with  the  second  condition.  The  simulation 
parameters  are  listed  below. 

The  start-up  and  transient  behaviors  are  analyzed  by  using 
these  conditions  and  the  simulated  results  are  discussed  in  the 
following  sections. 

The  operating  conditions  for  a  start-up  and  transient  response 
are  described  as  follows: 

•  The  start-up  temperature  for  the  cells  is  initially  set  to 
298.15  K.  A  look-up  table  was  created  for  the  lowest  mem¬ 
brane  temperature  among  the  10  cells  as  an  input  and  the 
current  density  as  an  output.  The  value  of  current  density  con¬ 
tinuously  follows  the  increase  of  the  temperature  in  the  mem¬ 
brane  in  order  to  quickly  raise  the  temperature  to  353 . 1 5  K  that 
is  assumed  as  a  typical  working  temperature.  In  addition,  a 
built-in  controller  regulates  the  temperature  in  the  membrane 
to  prevent  dehydration.  Whenever  the  temperature  in  differ¬ 
ent  cells  exceeds  the  desired  working  temperature,  a  coolant 
subsystem  is  turned  on  and  off  to  extract  the  excessive  heat 
produced  in  the  stack. 

•  A  load  current  with  multiple  steps  is  applied  to  the  stack  model 
and  the  transient  responses  analyzed. 

3.3.1.  Start-up 

3. 3. 1.1.  Case  1:  vapor  phase  ( low  cathode  inlet  RH  and  fast 
water  removal).  Fig.  3  shows  the  current  load,  the  transient 
behavior  of  the  temperature  for  the  membrane  layer  of  differ¬ 
ent  cells,  as  well  as  the  stack  temperature  profile  at  the  50,  360 
and  450  s  during  the  start-up.  In  the  given  start-up  condition 
(Table  4);  it  took  more  than  6  min  for  a  membrane  layer  with  the 
highest  temperature  to  reach  the  desired  working  temperature 
(Fig.  3b). 

Generally,  the  temperature  of  each  layer  rises  during  a  start¬ 
up.  Due  to  ohmic  losses  by  the  membrane  resistances  and  the 

Table  4 


working  condition  of  the  stack 


Gas  inlet  condition 

Pin 

(Pin 

RHin 

fo2  or  CO 

Cathode 

141000 

2.5 

0.8 

0.21 

Anode 

141000 

1 

0.001  or  0.5 

0 

heat  released  by  the  chemical  reaction  in  the  cathode  catalysts, 
the  membranes  and  cathode  catalysts  show  the  highest  value 
among  others.  In  contrast,  the  losses  on  the  anodic  side  are  neg¬ 
ligible.  Therefore,  the  high  losses  of  cathodic  catalysts  and  the 
corresponding  heat  generated  leads  to  an  asymmetrical  temper¬ 
ature  distribution  throughout  the  cell  [15,17,23]. 

In  fact,  the  temperature  of  membranes  rises  to  a  reference 
temperature  with  the  associated  rising  time  strongly  influenced 
by  the  location  of  a  cell  in  the  stack.  As  a  result,  the  total  heat 
capacity  of  an  endplate  assembly  is  much  higher  than  the  one 
of  the  layers  in  the  cells.  The  heat  generated  in  the  end  cells 
are  quickly  transferred  to  endplates  and  stored  there  rather  than 
stored  in  the  cell  itself.  The  temperature  of  the  membrane  layers 
in  the  middle  of  the  stack  increases  more  rapidly  than  the  end 
ones.  For  example,  the  5th  cell  shows  the  highest  temperature, 
while  the  10th  cell  shows  the  lowest  (Fig.  3b). 

The  profile  of  the  stack  temperature  in  Fig.  3c-e  shows  dif¬ 
ferent  dynamics  with  an  asymmetric  shape  that  differs  from  the 
one  with  symmetric  shape  [9]  or  smooth  curve  [20] .  In  addition, 
an  unsteady  and  abrupt  drop  of  temperature  at  the  end  cells  is 
observed,  caused  by  the  difference  between  ambient  and  cell 
temperature.  The  large  heat  capacity  of  the  bus  plate  and  the 
small  heat  transfer  coefficient  of  the  endplate  lead  to  a  negligi¬ 
ble  heat  transfer. 

At  the  50th  second  after  the  start-up  (Fig.  3c),  the  nine  cells 
show  the  temperature  walls  at  MEA  except  the  cell  number 
1.  In  fact,  the  cathode  catalyst  layer  of  the  cell  1  is  located 
near  the  endplate  and  possibly  conducts  a  large  amount  of 
heat  generated  toward  the  endplates.  Compared  to  the  first 
cell,  the  heat  generated  at  the  cathode  catalyst  in  the  10th 
cell  is  being  kept  by  two  more  layers  of  the  membrane  and 
anodic  catalysts.  Thus,  the  temperature  of  the  cathode  catalysts 
layer  in  the  first  cell  becomes  lower  than  the  one  of  the  10th 
cell. 

At  the  360th  second  (Fig.  3d),  the  membrane  of  the  fifth 
cell  shows  the  highest  temperature  at  80  °C,  while  six  cells  still 
show  a  temperature  wall  with  an  asymmetric  distribution.  The 
temperature  wall  is  generated  in  the  10th  cell  before  blocking 
the  conduction  of  the  heat  generated  by  the  rest  of  the  cells. 
The  transfer  of  heat  occurs  through  the  left  side  of  cells.  As  a 
result,  the  temperature  of  the  left  end  plate  assembly  rises  more 
than  the  one  of  the  right  end  plate  assembly.  Subsequently,  the 
temperature  walls  of  cells  located  at  the  right-hand  side  of  the 
stack  disappear.  Then,  the  temperature  of  the  cell  rises  gradually, 
while  the  first  cell  does  continuously.  Finally,  the  temperature  of 
the  cathode  catalysts  layer  in  the  10th  cell  becomes  lower  than 
the  one  of  the  first  cell. 

At  the  450th  second,  the  numbers  of  temperature  walls 
increase  again.  As  the  temperature  at  the  end  plate  assembly 
increases,  the  gradient  of  the  cells  and  the  end  plates  becomes 
lower  and  thus  the  heat  conduction  from  the  cells  decrease.  As  a 
result,  the  heat  generated  accumulates  in  the  end  cells,  while  the 
heat  in  the  inner  cells  causes  a  temperature  rise  in  the  coolant 
that  offsets  the  end  cells  with  a  temperature  rise.  Finally,  the 
temperature  of  the  stack  becomes  more  uniform.  The  middle 
cells  reach  a  steady  state,  while  the  temperature  dynamics  of  the 
end  cells  have  not  completed  yet. 
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Fig.  4  shows  a  dynamic  behavior  of  the  output  voltage  of 
the  individual  cell  and  output  voltage  referred  to  the  fifth  cell 
voltage  that  shows  the  highest  value. 

Generally,  the  cells  located  in  the  inner  side  of  a  stack  show 
low  losses  of  over-potentials  because  of  the  high  temperature. 
The  middle  cells  illustrate  a  higher  output  voltage  than  those  of 
the  end  cells. 

The  behavior  of  the  two  end  cells  is  particularly  different 
from  the  rest  of  the  cells  (Fig.  4a).  The  voltage  of  all  the  cells 
tends  to  follow  the  increase  of  temperature  but  those  of  end  cells 
show  a  decrease.  In  fact,  the  lower  temperature  at  the  end  cells 
slows  down  the  removal  of  water  continuously  generated,  which 
blocks  the  influx  of  the  reactants  toward  the  catalysts  through  the 
GDL.  Likewise,  water  generated  in  the  catalyst  is  hard  to  remove, 
which  causes  a  high  over-potential.  This  simulation  result  is 
comparable  to  the  experimental  results  [21]. 

The  output  voltage  of  the  cells  is  shown  in  Fig.  4b.  The  volt¬ 
age  increases  until  300  s  and  remains  almost  constant  affected 
mainly  by  dehydrated  membranes  and  the  associated  low  proton 
conductivities.  On  the  other  hand,  the  10th  cell  shows  less  loss  in 
the  over-potential  than  the  first  cell  until  300  s,  because  the  tem¬ 
perature  of  the  cathode  catalyst  at  the  first  cell  becomes  larger 
than  the  one  at  the  10th  cell.  As  soon  as  the  coolants  start  control¬ 
ling  temperature  around  350  s,  the  output  voltage  of  cells  except 
two  end  cells  starts  decreasing  until  the  current  reaches  around 
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Fig.  5.  Dynamic  behavior  of  membrane  water  uptake  for  a  stack. 


at  the  right-hand  side  show  a  higher  rate  of  water  uptake  than  its 
counterpart  in  the  stack.  The  differences  in  temperature  between 
the  cathodic  electrodes  of  cells  at  the  left-hand  side  of  the  stack 
and  at  the  right-hand  side  leads  to  different  saturating  pressures 
(Eq.  (9)),  although  the  amount  of  water  generated  is  the  same 
at  the  load  current.  The  lower  the  temperature  of  the  cathode 
gas,  the  lower  the  saturation  pressure.  Thus,  the  outlet  gas  car¬ 
ries  less  water  vapor  out  of  cells,  and  consequently,  the  water 
concentration  in  the  membrane  becomes  higher. 


6.853e  -  1Ta  -  7.432e  -  4 T3  +  0.304r2  -  55.613T  +  3831.8  if  (T  <  333.15  K) 

6.853e  -  1Ta  -  7.432e  -  4 T3  +  0.304r2  -  55.613T  +  3831.8  if  (333.15  K  <  T  <  433.15  K) 


0.6  A  cm-2.  The  operating  condition  of  the  anodic  inlet  RH  is 
intentionally  set  to  0.6  from  0.001  in  order  to  study  the  effects  of 
the  humidity  on  the  membranes.  As  a  result,  the  voltage  of  the 
cells  starts  increasing  because  of  the  increased  membrane  proton 
conductivity.  The  recovery  speed  of  membrane  conductivity  of 
the  individual  cell  becomes  dependent  on  the  cell  location  in 
the  stack,  whose  physical  reasons  are  described  in  detail  in  the 
following  section. 

Fig.  4b  shows  percentage  rates  of  the  cell  voltages.  The  volt¬ 
age  of  the  end  cells  amounts  to  half  of  the  voltage  produced  by 
the  central  cells.  The  derivation  of  the  central  cells  varies  at  a 
rate  of  around  20%,  which  also  worsens  overall  performance  of 
the  stack. 

The  proton  conductivity  of  membrane  currently  available 
strongly  depends  upon  water  content.  Therefore,  an  analysis  is 
performed  to  study  the  effect  of  the  anodic  RH  and  cell  loca¬ 
tion  on  membranes.  The  inlet  RH  has  been  drastically  changed 
from  0.001  to  0.6  at  the  380  s.  Fig.  5  shows  simulated  results  of 
the  water  uptake  of  membrane  layers  in  the  10  cells  during  the 
start-up. 

At  the  low  anode  inlet  gas  humidification,  water  present  in  the 
membrane  is  only  influenced  by  water  generated  in  the  catalysts. 
Thus,  the  water  uptake  continuously  decreases  and  correspond¬ 
ing  proton  conductivity  gets  lower. 

At  an  application  of  a  high  humidification,  the  water  uptake 
in  the  membranes  of  all  the  cells  increases.  Particularly,  the  cells 


A  comparison  of  water  uptake  shows  a  large  difference 
affected  by  the  temperature  distribution  in  the  stack.  For  exam¬ 
ple,  the  central  cells  become  quickly  hydrated  in  comparison  to 
the  end  cells.  Under  a  fully  consumed  fuel,  the  anode  gas  gets 
saturated.  Therefore,  the  electro-osmotic  influence  at  the  anodic 
boundary  of  each  cell  on  their  membranes  becomes  the  same. 
However,  the  temperature  distribution  leads  to  much  higher 
water  concentration  and  diffusion  coefficient  in  membrane  lay¬ 
ers  at  the  central  cells  than  the  ones  at  the  end  cells. 

If  the  influence  of  the  liquid  water  is  neglected,  the  diffusion 
of  water  from  the  anodic  electrodes  to  membranes  in  the  central 
cells  becomes  higher  than  the  one  in  the  end  cells.  Conversely, 
the  humidity  at  the  cathodic  inlet  is  low  and  the  temperature  is 
high,  so  the  gas  is  not  saturated.  Due  to  the  temperature  differ¬ 
ence  in  the  cells,  the  saturation  pressure  of  end  cells  is  lower 
than  the  one  of  the  central  cells  (Eq.  (9))  Consequently,  water 
activity  at  the  cathodic  sides  of  the  end  cells  is  higher,  which 
makes  it  easier  to  transport  water  out  of  membranes  by  electro- 
osmotic  force  at  the  cathode  boundary.  As  a  result,  the  recovery 
speed  of  water  contents  in  the  central  cells  is  generally  faster 
than  the  one  in  the  end  cells.  However,  the  low  water  concen¬ 
tration  in  the  cathode  sides  of  the  central  cells  accelerates  water 
diffusion  from  membranes  to  cathodes,  which  again  slows  down 
the  recovery  speed. 

3. 3. 1.2.  Liquid  phase  (high  cathode  inlet  RH  and  slow  water 
removal ).  Fig.  6  shows  the  current  load,  the  transient  behavior  of 
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Fig.  8.  Dynamic  membrane  water  uptake  of  the  10  cells. 


the  temperature  for  the  membrane  layers  of  different  cells  as  well 
as  the  stack  temperature  profile  at  the  50,  270  and  350  s.  Even 
with  the  same  lookup  table  of  current  and  membrane  tempera¬ 
ture  for  the  start-up,  the  speed  of  the  start-up  becomes  quicker 
than  the  one  in  the  first  case.  It  takes  about  4.5  min  for  the  mem¬ 
brane  layer  with  the  highest  temperature  among  others  to  reach 
the  desired  working  temperature.  In  fact,  the  water  generated 
in  the  cell  is  removed  as  liquid  phase,  so  the  latent  heat  neces¬ 
sary  for  getting  water  evaporated  is  stored  in  the  stack.  Fig.  6b 
shows  the  temperature  waveforms  of  membranes.  A  comparison 
between  Figs.  3b  and  6b  shows  two  differences:  (1)  tempera¬ 
ture  of  membranes  continuously  rises,  even  if  the  coolant  circuit 
turns  on  and  (2)  the  gap  between  the  membrane  temperature 
of  the  5th  cell  and  the  10th  cell  is  calculated  by  the  time  the 
coolant  circuit  turns  on  75  s  later.  The  difference  in  both  gaps 
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has  been  decreased  from  8  °C  in  the  first  case  to  2  °C  in  this 
case. 

Therefore,  the  low  RH  at  cathode  inlet  in  the  first  case  enables 
the  coolant  to  easily  control  the  stack  temperature  at  the  start-up. 

Fig.  6  d-f  show  the  stack  temperature  profile  at  the  50, 270  and 
350  s.  A  comparison  with  the  first  case  shows  relatively  higher 
walls  of  temperature  with  increased  numbers  in  the  stack. 

Fig.  7  shows  the  dynamic  response  of  the  output  voltage  of 
the  10  cells  on  the  load  current  and  the  corresponding  percent¬ 
age  rate.  A  comparison  with  the  first  case  shows  a  low  output 
voltage  at  about  0.6  A  cm-2,  even  though  the  temperature  of  the 
stack  and  membrane  conductivity  are  higher  than  before.  As  a 
matter  of  fact,  the  liquid  water  generated  at  the  catalysts  layers 
blocks  the  influx  of  oxygen  and  finally  leads  to  a  starvation  of 
oxygen. 

Fig.  8  shows  the  water  uptake  in  the  membrane  layer  of  the 
10  cells  at  the  start-up.  The  humidified  cathode  inlet  gas  slightly 
slows  down  dehydration  of  the  membranes  until  the  300  s  when 
the  anode  inlet  gas  starts  to  be  humidified  from  0.00 1  to  0.6.  Then 
the  water  uptake  of  the  membrane  stops  decreasing.  A  compari¬ 
son  with  the  results  from  the  changed  cathode  and  anode  inlet  gas 
RH  shows  that  water  transport  at  the  boundary  of  membranes  is 
dominated  by  the  electro-osmotic  drag  force.  Therefore,  humid¬ 
ification  on  the  cathode  inlet  gas  is  not  sufficient  to  prevent 
dehydration,  which  requires  additional  humidification  of  the  gas 
at  the  anodic  inlet. 


3.3.2.  Transient  response 

3. 3. 2.1.  Case  1:  vapor  phase  ( low  cathode  inlet  RH  and  fast 
water  removal).  Fig.  9  shows  the  load  current  and  the  simulated 
results  for  the  temperature  of  the  membrane  layers  and  the  stack 
at  the  580,  650  and  800  s. 

When  the  current  drops  with  0.15  and  0.3  A  cm-2,  the  tem¬ 
perature  of  all  membranes  drops  by  1°  and  1.5°,  respectively 
(Fig.  9b).  Primarily,  the  heat  associated  with  the  current  density 
is  responsible  for  the  rise  or  drop  in  temperature.  Therefore,  the 
recovery  behavior  and  final  value  is  strongly  influenced  by  the 
amount  of  heat  dependent  upon  the  amplitude  of  current  density. 
For  example,  a  step  with  a  high  current  density  (0.45  A  cm-2) 
causes  an  immediate  rise  in  temperature,  while  a  low  current 
density  (0.15  A  cm-2)  takes  100  s  for  the  central  membrane  to 
reach  a  desired  working  temperature  of  353.15  K. 

The  simulation  shows  that  the  temperature  rise  of  the  end  cells 
at  the  high  current  density  is  faster  than  the  one  at  the  low  current. 
Besides,  the  heat  depends  upon  the  current  density  because  the 
coolant  controlling  for  a  working  temperature  carries  over  the 
heat  from  the  central  cells  and  offsets  the  temperature  of  the 
end  cells.  Similarly,  the  transfer  of  the  heat  to  the  end  cells 
gets  smaller.  When  the  amplitude  of  the  current  is  small,  the 
corresponding  heat  generated  in  the  central  cell  becomes  smaller. 
Those  effects  are  illustrated  in  Figs.  9d-f. 

Fig.  10  shows  the  dynamic  response  of  the  cell  output  volt¬ 
age  and  membrane  water  uptake  of  the  10  cells  at  the  same  load 


450  500  500  600  650  700  750  800  850  900 


(a)  Time  (sec) 


a> 

CD 

CO 


CD 

O 

i_ 

0) 

a. 


100 


80 


■  Cell  1 
o  Cell  2 
^  Cell  3 
v  Cell  4 

*  Cell  5 
a  Cell  6 

*  Cell  7 

*  Cell  8 

*  Cell  9 
«  Cell  10 


500  600  700  800  900 


Time  (sec) 


(d) 


Fig.  10.  (a  and  b)  Dynamic  voltage  output  of  the  10  cells,  (c)  Voltage  output  percentage  of  the  10  cells,  (d)  Dynamic  membrane  water  uptake  of  the  10  cells. 


284 


Y.  Shan,  S.-Y.  Choe  /  Journal  of  Power  Sources  158  (2006)  274-286 


current.  The  overshoot  of  voltage  is  caused  by  several  physical 
factors:  (1)  proton  concentration  at  the  cathode  catalysts  [16], 
(2)  oxygen  concentration  at  the  catalyst  [22],  and  (3)  variation  of 
the  cell  working  temperature.  Temperature  influences  the  water 
saturation  pressure  (Eq.  (9))  and  the  effective  areas  of  cathode 
catalysts  as  well  as  oxygen  concentration.  More  details  are  dis¬ 
cussed  in  Section  3. 3. 2. 2. 

Fig.  10b  and  c  shows  responses  of  water  uptake  on  the  load 
current.  The  operating  condition  for  the  RH  at  the  anode  inlet  is 
set  to  0.001  from  0.6  at  about  550  s  when  the  current  density  is 
lower  than  0.6  A  cm-2.  The  dependency  of  the  membrane  water 
uptake  on  the  temperature  is  identical  with  the  previous  analyses. 

33.2.2.  Case  2:  liquid  phase  (high  cathode  inlet  RH  and  slow 
water  removal ).  The  step  load  current  has  been  applied  to  study 


the  effects  of  liquid  water  removal  on  the  dynamics  (Fig.  11a). 
Fig.  1  lb-f  show  the  dynamic  temperature  behavior  of  the  mem¬ 
brane  layers  and  the  stack  temperature  profile  at  the  400,  550 
and  800  s.  A  comparison  with  the  first  case  of  operating  condi¬ 
tion  shows  a  high  temperature  drop  in  the  stack.  It  is  caused  by 
the  large  amount  of  heat  generated  in  the  catalysts  when  water 
is  generated  and  removed  as  liquid.  Consequently,  the  temper¬ 
ature  cannot  be  controlled  with  the  set  coolant  flow  rate  until 
the  current  density  is  decreased.  At  round  900  s,  the  temperature 
difference  between  the  membrane  of  the  central  cell  and  two 
end  cell  amounts  to  about  5  and  13  °C,  respectively,  less  than 
the  previously  calculated  values  of  5.5  and  15  °C  even  at  the 
same  amplitude  of  the  current  load. 

Generally,  the  amplitude  of  the  peak  temperature  at  both 
high  and  low  current  density  is  higher  than  the  first  case  and 
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Fig.  11.  (a)  Step  current  input,  (b)  Dynamic  membrane  temperature  of  the  10  cells,  (c)  Stack  temperature  distribution  at  the  400  s.  (d)  Stack  temperature  distribution 
at  the  550  s.  (e)  Stack  temperature  distribution  at  the  800s. 
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Fig.  12.  (a  and  b)  Dynamic  voltage  output  of  the  10  cells,  (c)  Voltage  output  percentage  of  the  10  cells,  (d)  Dynamic  membrane  water  uptake  of  the  10  cells. 


the  numbers  of  temperature  walls  at  MEA  is  more  than  before 
(Fig.  lld-f). 

Fig.  12  shows  the  same  dynamic  responses  of  output  voltage 
of  all  the  cells.  The  amplitude  of  the  overshoot  of  the  output  cell 
voltage  is  much  larger  than  the  first  case  according  to  Fig.  7a, 
which  is  strongly  dependent  on  the  temperature.  An  analysis 
shows  that  a  sudden  temperature  drop  with  a  high  current  causes 
high  amplitude  of  the  overshoot.  Conversely,  at  the  low  current 
step,  the  heat  generated  can  be  more  effectively  extracted  by 
the  same  condition  of  the  coolant  circuit,  so  the  influence  of 
temperature  on  the  voltage  overshoot  becomes  less. 

In  addition,  the  waveform  of  voltage  of  the  10  cells  does  not 
have  the  same  tendency  but  fluctuates  with  distortions  at  special 
operating  conditions.  For  example,  the  output  voltage  of  the  sec¬ 
ond  cell  at  an  instant  can  be  higher  than  the  fifth  cell  expected 
with  the  highest  voltage.  The  open  circuit  voltage  (OCV)  is  a 
function  of  the  temperature  whose  derivative  shows  a  negative 
value.  Therefore,  OCV  decreases  when  the  temperature  rises. 
Consequently,  the  second  OCV  is  higher  than  the  fifth.  In  addi¬ 
tion,  the  high  temperature  causes  more  water  to  be  stored  as 
vapor  at  the  electrodes.  Particularly,  the  second  cell  has  lower 
temperature  than  the  one  in  the  fifth,  which  leads  to  a  less  water 
vapor.  The  concentration  of  oxygen  becomes  larger  and  the  asso¬ 
ciated  over-potential  does  smaller. 

On  the  other  hand,  temperature  increase  accelerates  dehydra¬ 
tion  of  the  membranes  (Fig.  12c)  that  decreases  the  membrane 
proton  conductivity.  Therefore,  corresponding  proton  conduc¬ 


tivity  of  the  second  cell  is  higher  than  the  fifth  cell.  Although 
the  liquid  water  in  the  fifth  cell  is  less  than  the  one  in  the  second, 
the  overall  cell  performance  of  the  second  cell  is  better  because 
of  the  three  reasons  aforementioned. 

4.  Conclusion 

This  paper  presents  a  model  for  the  PEM  fuel  cell  stack  based 
on  single  cells  with  layers.  Simulation  is  conveyed  to  analyze 
the  effects  of  temperature  on  static  and  dynamic  characteris¬ 
tics  of  cells  and  finally  the  stack.  Two  operating  conditions 
are  considered  to  investigate  behaviors  of  the  stack  at  a  start¬ 
up  and  a  step  load.  The  results  show  asymmetrical  temperature 
distributions  through  the  stack,  varying  dynamically  with  oper¬ 
ating  conditions  and  applied  loads.  The  analyses  suggest  a  new 
consideration  for  design  and  operating  criteria  to  improve  the 
performance  of  cells  and  optimize  components  of  the  power 
system.  In  particular,  proper  thermal  management  is  required  to 
prevent  dehydration  of  membranes  and  secure  longevity  of  the 
cells.  Fikewise,  optimal  operations  of  the  stack  can  be  derived 
by  this  simulation  and  associated  new  control  strategies  can  be 
developed.  Major  findings  are  summarized  as  follows: 

•  Temperature  distribution  through  the  cells  is  asymmetric 
because  of  the  heat  generated  in  the  cathode  sides  are  higher 
than  in  the  anodic  sides.  This  asymmetric  effect  coupled  with 
the  end  plate  finally  leads  to  an  asymmetrical  temperature 
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profile  in  the  stack  (Figs.  3b,  6b,  9b  and  lib)  determined 
by  thermal  conductivities  of  the  end  plates.  The  temperature 
wall  in  the  cells  blocks  further  conduction  of  heat  generated 
by  other  cells  contingent  on  their  location,  and  strongly  influ¬ 
ences  the  static  and  dynamic  characteristics; 

•  The  heat  conductivity  of  the  membrane  varies  during  oper¬ 
ations  because  of  the  continuous  variation  of  water  content 
and  the  consequent  swelling  phenomena,  which  results  in  an 
increase  of  electrical  and  thermal  resistance. 

•  The  latent  heat  of  water  produced  at  the  catalysts  by  the  reac¬ 
tions  can  be  stored  in  the  stack  and  requires  a  proper  cooling 
strategy  at  a  start-up  when  the  cathode  outlet  gas  does  not 
sufficiently  carry  all  the  water  vapor.  In  this  case,  the  time 
necessary  for  a  start-up  is  relatively  short.  Thus,  the  temper¬ 
ature  difference  between  the  central  cells  and  the  end  cells  is 
larger  than  the  one  with  water  removal  during  the  vapor  phase. 

•  Non-uniform  temperature  distribution  through  the  stack  can 
be  minimized  by  properly  coupling  the  coolant  for  the  central 
cells  with  the  end  cells  which  can  offset  the  temperature  of 
the  end  cells. 

5.  Future  work 

We  plan  to  add  an  accurate  mass  transport  model  to  the  current 
stack  model  to  study  the  dynamic  effects  of  fuels  and  reactants  on 
the  stack.  The  models  will  then  be  integrated  into  the  exiting  sys¬ 
tem  model  and  implemented  on  a  real-time  system  to  optimize 
components  of  a  fuel  cell  power  system  and  design  advanced 
controls.  Further  improvements  are  anticipated  by  taking  into 
account  the  gas  dynamics  in  the  flow  channel.  At  the  same  time, 
we  plan  to  validate  the  models  in  close  collaboration  with  indus¬ 
tries  and  are  looking  for  a  partner  to  work  with. 
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